clear all


use "$root\data\use.dta", replace


capture erase "$root/results_replication/omit/ovb_re.xls"
capture erase "$root/results_replication/omit/ovb_re.txt"
capture erase "$root/results_replication/omit/ovb_re_corn.xls"
capture erase "$root/results_replication/omit/ovb_corn_re.txt"
capture erase "$root/results_replication/omit/ovb_soybeans_re.xls"
capture erase "$root/results_replication/omit/ovb_soybeans_re.txt"

**# Omitted variables with drought and wetness
*First get results with wetness constrols from the main result and ones with omitted variables

foreach crop in corn soybeans {
	local c_name = strproper("`crop'")

	preserve
	local c = cond("`crop'" == "corn", "c", "s")
	rename `c'_d*_* d*_*
	rename `c'_tot_* tot_*
	rename `c'_pdsi_w_* pdsi_w_*
	rename `c'_pmdi_w_* pmdi_w_*
	rename `c'_tbin*_* tbin*_*
	rename `c'_prec_* prec_*
	rename `crop'_production00 production00

**## Harvest ratio
	reghdfe `crop'_hratio d*_b d*_p d*_g pdsi_w_b tbin*b pdsi_w_p tbin*p pdsi_w_g tbin*g if irrigated_ratio<0.1 & N_year_`crop'>10 & lng>-100, absorb(fips i.state#c.year i.state#c.year2 i.year) vce(cluster fips stateyear)
	estimate store `crop'_hratio_wet
	
	reghdfe `crop'_hratio d*_g pdsi_w_g tbin*g if irrigated_ratio<0.1 & N_year_`crop'>10 & lng>-100, absorb(fips i.state#c.year i.state#c.year2 i.year) vce(cluster fips stateyear)
	estimate store `crop'_hratio_omit

**## Yield
	reghdfe ln`crop'_yield d*_b d*_p d*_g pdsi_w_b tbin*b pdsi_w_p tbin*p pdsi_w_g tbin*g if irrigated_ratio<0.1 & N_year_`crop'>10 & lng>-100, absorb(fips i.state#c.year i.state#c.year2 i.year) vce(cluster fips stateyear)
	estimate store `crop'_yield_wet
	
	reghdfe ln`crop'_yield d*_g pdsi_w_g tbin*g if irrigated_ratio<0.1 & N_year_`crop'>10 & lng>-100, absorb(fips i.state#c.year i.state#c.year2 i.year) vce(cluster fips stateyear)
	estimate store `crop'_yield_omit
		
**## Production
	reghdfe ln`crop'_production d*_b d*_p d*_g pdsi_w_b tbin*b pdsi_w_p tbin*p pdsi_w_g tbin*g if irrigated_ratio<0.1 & N_year_`crop'>10 & lng>-100, absorb(fips i.state#c.year i.state#c.year2 i.year) vce(cluster fips stateyear)
	estimate store `crop'_production_wet
	
	reghdfe ln`crop'_production d*_g pdsi_w_g tbin*g if irrigated_ratio<0.1 & N_year_`crop'>10 & lng>-100, absorb(fips i.state#c.year i.state#c.year2 i.year) vce(cluster fips stateyear)
	estimate store `crop'_production_omit

restore
}


**Results with misspecification: using whole year drought
foreach crop in corn soybeans {
	local c_name = strproper("`crop'")

	preserve
	local c = cond("`crop'" == "corn", "c", "s")
	rename `c'_d*_* d*_*
	rename `c'_tot_* tot_*
	rename `c'_pdsi_w_* pdsi_w_*
	rename `c'_pmdi_w_* pmdi_w_*
	rename `c'_tbin*_* tbin*_*
	rename `c'_prec_* prec_*
	rename `crop'_production00 production00
	
	forvalues s = 0/3 {
		replace d`s'_g = (d`s'_b_w + d`s'_p_w + d`s'_g_w)*7/365 //What is this for? Whole year drought?
	}

	reghdfe `crop'_hratio d*_g pdsi_w_g tbin*g if irrigated_ratio<0.1 & N_year_`crop'>10 & lng>-100, absorb(fips i.state#c.year i.state#c.year2 i.year) vce(cluster fips stateyear)
	estimate store `crop'_hratio_mis
	
	reghdfe ln`crop'_yield d*_g pdsi_w_g tbin*g if irrigated_ratio<0.1 & N_year_`crop'>10 & lng>-100, absorb(fips i.state#c.year i.state#c.year2 i.year) vce(cluster fips stateyear)
	estimate store `crop'_yield_mis
	
	reghdfe ln`crop'_production d*_g pdsi_w_g tbin*g if irrigated_ratio<0.1 & N_year_`crop'>10 & lng>-100, absorb(fips i.state#c.year i.state#c.year2 i.year) vce(cluster fips stateyear)
	estimate store `crop'_production_mis
	
	restore
}

foreach outcome in hratio yield production {
	local y_name = cond("`outcome'" == "hratio", "Harvested Ratio", strproper("`outcome'"))
	local order = cond("`outcome'" == "hratio", "a", cond("`outcome'" == "yield", "b", "c"))
	coefplot  /// 
		(corn_`outcome'_omit, label(Growing Season Only) offset(-0.1) ciopts(color(black) lpattern(shortdash))  msymbol(T)) ///
		(corn_`outcome'_wet, label(All Seasons Separately) offset(0) ciopts(color(black)) msymbol(S)) ///
		(corn_`outcome'_mis, label(Whole Year) offset(0.15) ciopts(color(black) lpattern(dot))), bylabel(Corn) || /// 
		(soybeans_`outcome'_omit, label(Growing Season Only) offset(-0.1) ciopts(color(black) lpattern(shortdash)) msymbol(T)) /// 
		(soybeans_`outcome'_wet, label(All Seasons Separately) offset(0) ciopts(color(black)) msymbol(S)) ///
		(soybeans_`outcome'_mis, label(Whole Year) offset(0.15) ciopts(color(black) lpattern(dot))), bylabel(Soybeans) ||, ///
		mcolor(black) mfcolor(white) vertical keep(d*_g) yline(0, lcolor(gray) lwidth(0.2)) /// 
		byopts(compact cols(1) graphregion(fcolor(white))) subtitle(, fcolor(white) bmargin(top) bcolor(white)) ///
		groups(d*_g = "Drought Severity") ytitle("ln(`y_name')") ylab(,labs(small)) xlabel(1 "D0" 2 "D1" 3 "D2" 4 "D3") ///
		legend(rows(1))

}


foreach crop in corn soybeans {
	foreach outcome in hratio yield production {

	local c_name = strproper("`crop'")

	preserve
	local c = cond("`crop'" == "corn", "c", "s")
	rename `c'_d*_* d*_*
	rename `c'_tot_* tot_*
	rename `c'_pdsi_w_* pdsi_w_*
	rename `c'_pmdi_w_* pmdi_w_*
	rename `c'_tbin*_* tbin*_*
	rename `c'_prec_* prec_*
	rename `crop'_production00 production00

	if "`outcome'" == "hratio" {
		qui reghdfe `crop'_`outcome' tot_g tbin*_g pdsi_w_g if irrigated_ratio<0.1 & N_year_`crop'>10 & lng>-100, absorb(fips i.state#c.year i.state#c.year2 i.year) vce(cluster fips stateyear)
	}
	else {
		qui reghdfe ln`crop'_`outcome' tot_g tbin*_g pdsi_w_g if irrigated_ratio<0.1 & N_year_`crop'>10 & lng>-100, absorb(fips i.state#c.year i.state#c.year2 i.year) vce(cluster fips stateyear)
	}
	
	if "`outcome'" == "hratio" {
		qui outreg2 using "$root/results_replication/tables/ovb_`crop'_re.xls", excel bdec(4) sdec(4) keep(tot_g) nocon adjr2 replace
	}
	else {
		qui outreg2 using "$root/results_replication/tables/ovb_`crop'_re.xls", excel bdec(4) sdec(4) keep(tot_g) nocon adjr2 
	}

	if "`outcome'" == "hratio" {
		qui reghdfe `crop'_`outcome' tot_b tot_p tot_g tbin*_b pdsi_w_b tbin*_p pdsi_w_p tbin*_g pdsi_w_g if irrigated_ratio<0.1 & N_year_`crop'>10 & lng>-100, absorb(fips i.state#c.year i.state#c.year2 i.year) vce(cluster fips stateyear)
	}
	else {
		qui reghdfe ln`crop'_`outcome' tot_b tot_p tot_g tbin*_b pdsi_w_b tbin*_p pdsi_w_p tbin*_g pdsi_w_g if irrigated_ratio<0.1 & N_year_`crop'>10 & lng>-100, absorb(fips i.state#c.year i.state#c.year2 i.year) vce(cluster fips stateyear)
	}
	qui outreg2 using "$root/results_replication/tables/ovb_`crop'_re.xls", excel bdec(4) sdec(4) keep(tot_g tot_b tot_p) nocon adjr2

	* Equality tests
	if "`outcome'" == "hratio" {
		qui reghdfe `crop'_`outcome' tot_g tbin*_g pdsi_w_g if irrigated_ratio<0.1 & N_year_`crop'>10 & lng>-100, absorb(fips i.state#c.year i.state#c.year2 i.year) cluster(fips stateyear)
		qui keep if e(sample)
		qui expand 2, gen(sample)
		qui gen cons = 1
		qui reghdfe `crop'_`outcome' i.sample#c.(tot_g tbin*_g pdsi_w_g cons) c.sample#c.(tot_b tbin*_b pdsi_w_b tot_p tbin*_p pdsi_w_p ) if irrigated_ratio<0.1 & N_year_`crop'>10 & lng>-100, absorb(sample#fips sample#i.state#c.year sample#i.state#c.year2 sample#i.year) cluster(fips stateyear)
		}
	else {
		qui reghdfe ln`crop'_`outcome' tot_g tbin*_g pdsi_w_g if irrigated_ratio<0.1 & N_year_`crop'>10 & lng>-100, absorb(fips i.state#c.year i.state#c.year2 i.year) cluster(fips stateyear)
		qui keep if e(sample)
		qui expand 2, gen(sample)
		qui gen cons = 1
		qui reghdfe ln`crop'_`outcome' i.sample#c.(tot_g tbin*_g pdsi_w_g cons) c.sample#c.(tot_b tbin*_b pdsi_w_b tot_p tbin*_p pdsi_w_p ) if irrigated_ratio<0.1 & N_year_`crop'>10 & lng>-100, absorb(sample#fips sample#i.state#c.year sample#i.state#c.year2 sample#i.year) cluster(fips stateyear)		
	}

	qui test 0.sample#c.tot_g=1.sample#c.tot_g
	di "F-value of equality test between tot_g for `crop'_`outcome':" r(F)
	di "p-value of equality test between tot_g for `crop'_`outcome':" r(p)
	
	if "`crop'" == "corn" & "`outcome'" == "hratio" {
		qui outreg2 using "$root/results_replication/tables/ovb_re.xls", excel bdec(4) sdec(4) keep(sample#c.tot_g) nocon noobs addstat(F-Stat, r(F), P-Value, r(p)) replace
	}
	else {
		qui outreg2 using "$root/results_replication/tables/ovb_re.xls", excel bdec(4) sdec(4) keep(sample#c.tot_g) nocon noobs addstat(F-Stat, r(F), P-Value, r(p))
	}

	restore		
	}
}
